function [ctrl, aux] = GuidanceOutLoop(t, tf, statec, statef, ctrlp, auxp)
global Vs d2r m S R0 Vup Vdown alphamax dAziUp dAziDown alphaLD ef
rad  = statec(1);
lon = statec(2);
lat = statec(3);
V = statec(4);
azi = statec(6);
radf  = statef(1);
lonf = statef(2);
latf = statef(3);
Vf = statef(4);
fpa = statec(5);
alphap = ctrlp(1);
absBankp = ctrlp(2);
signOfBankp = ctrlp(3);
absBank = auxp(1);
absBankf = auxp(2);
e = 1/rad - V^2/2;
kh = -0.04 + 2e-5*V*Vs;
eps1 = 1E-3;
eps2 = 1E-3;
%% AOA
if V > Vup
    alphac = alphamax;
elseif V > Vdown
    alphac = (alphamax - alphaLD)*(V - Vdown)/(Vup - Vdown) + alphaLD;
else
    alphac = alphaLD;
end
ctrl(1) = alphac;
%% magnitude of bank
% StogoExp = sqrt((rad-radf)^2 + (acos(sin(lat)*sin(latf) + cos(lat)*cos(latf)*cos(lonf-lon)))^2);
StogoExp = acos(sin(lat)*sin(latf) + cos(lat)*cos(latf)*cos(lonf-lon));
tgoExp = tf - t;
% 牛顿迭代法
[bankbase, absBankf, F, G, iter] = newtonIteration2(@(absBank0,absBankf)PredIntegral(e, absBank0, ...
    absBankf, statec, statef, signOfBankp, StogoExp, tgoExp), absBank, absBankf, 10, eps1, eps2);
% QEGC补偿
alt = (rad-1)*R0;
[~, sonic, ~, rho] = atmoscoesa(alt);
Ma = V*Vs/sonic;
[CL, ~] = getCLCD(alphaLD, Ma);
[~, CD] = getCLCD(alphap, Ma);
q = 0.5*rho*(V*Vs)^2;
L = q*S*CL/m;
D = q*S*CD/m;
dhc = - D * 2 * 7110 / (V*Vs);
dh = (V*Vs) * sin(fpa);
ctrl(2) = acos((L*cos(bankbase) + kh * (dhc - dh))/L);
%% sign of bank
aziLos = atan(sin(lonf-lon)/(cos(lat)*tan(latf) - sin(lat)*cos(lonf-lon)));
dAzi = azi - aziLos;
if dAzi > dAziUp
    signOfBank = -1;
elseif dAzi > dAziDown
    signOfBank = signOfBankp;
else
    signOfBank = 1;
end
ctrl(3) = signOfBank;
aux(1) = bankbase;
aux(2) = absBankf;
aux(3) = iter;
end